function lab_img = rgb2lab(rgb_img)
% Convert RGB 2 LAB
% TuanND
% 03/30
%Pre-Step 1: convert sRGB -> [0 1]
nrgb_img = normalized(rgb_img);
% Step 1: convert nomalized sRGB -> linear RGB
lrgb_img = nrgb2lrgb(nrgb_img);
% Step 2
xyz_img = lrgb2xyz(lrgb_img);
% Step 3
lab_img = xyz2lab(xyz_img);
end
%--------------------------------------------------------------------------
function lab_img = xyz2lab(xyz_img)
[rows cols channel] = size(xyz_img);
x_channel = xyz_img(:,:,1);
y_channel = xyz_img(:,:,2);
z_channel = xyz_img(:,:,3);

Xw = 0.950;
Yw = 1;
Zw = 1.089;

X = x_channel/Xw;
Y = y_channel/Yw;
Z = z_channel/Zw;

f_Xw(X > 0.008856) = X(X > 0.008856) .^ (1/3);
f_Xw(X <= 0.008856) =  7.787 * X(X <= 0.008856) + 16/116;

f_Yw(Y > 0.008856) = Y(Y > 0.008856) .^ (1/3);
f_Yw(Y <= 0.008856) =  7.787 * Y(Y <= 0.008856) + 16/116;

f_Zw(Z > 0.008856) = Z(Z > 0.008856) .^ (1/3);
f_Zw(Z <= 0.008856) =  7.787 * Z(Z <= 0.008856) + 16/116;

l_channel = 116 * f_Yw - 16;
a_channel = 500 * (f_Xw - f_Yw);
b_channel = 200 * (f_Yw - f_Zw);
l_channel = reshape(l_channel, rows, cols);
a_channel = reshape(a_channel, rows, cols);
b_channel = reshape(b_channel, rows, cols);
lab_img(:,:,1) = l_channel;
lab_img(:,:,2) = a_channel;
lab_img(:,:,3) = b_channel;
end
%--------------------------------------------------------------------------
function xyz_img = lrgb2xyz(lrgb_img)
[rows cols channel] = size(lrgb_img);
convert_matrix = [0.412 0.358 0.180;0.213 0.715 0.072; 0.019 0.119 0.951];
r_channel = lrgb_img(:,:,1);
g_channel = lrgb_img(:,:,2);
b_channel = lrgb_img(:,:,3);

r_channel = reshape(r_channel, 1, []);
g_channel = reshape(g_channel, 1, []);
b_channel = reshape(b_channel, 1, []);

temp_xyz_img = convert_matrix * [r_channel; g_channel; b_channel];
x_channel = reshape(temp_xyz_img(1,:), rows, cols);
y_channel = reshape(temp_xyz_img(2,:), rows, cols);
z_channel = reshape(temp_xyz_img(3,:), rows, cols);
xyz_img(:,:,1) = x_channel;
xyz_img(:,:,2) = y_channel;
xyz_img(:,:,3) = z_channel;
end
%--------------------------------------------------------------------------
function lrgb_img = nrgb2lrgb(nrgb_img)
[rows cols channel] = size(nrgb_img);
nr_channel = nrgb_img(:,:,1);
ng_channel = nrgb_img(:,:,2);
nb_channel = nrgb_img(:,:,3);

sr_channel(nr_channel <= 0.04045) = nr_channel(nr_channel <= 0.04045) / 12.92;
sr_channel(nr_channel > 0.04045) = ((nr_channel(nr_channel > 0.04045) + 0.055) / 1.055).^2.4;
sr_channel = reshape(sr_channel, rows, cols);
sg_channel(ng_channel <= 0.04045) = ng_channel(ng_channel <= 0.04045) / 12.92;
sg_channel(ng_channel > 0.04045) = ((ng_channel(ng_channel > 0.04045) + 0.055) / 1.055).^2.4;
sg_channel = reshape(sg_channel, rows, cols);
sb_channel(nb_channel <= 0.04045) = nb_channel(nb_channel <= 0.04045) / 12.92;
sb_channel(nb_channel > 0.04045) = ((nb_channel(nb_channel > 0.04045) + 0.055) / 1.055).^2.4;
sb_channel = reshape(sb_channel, rows, cols);
lrgb_img(:,:,1) = sr_channel;
lrgb_img(:,:,2) = sg_channel;
lrgb_img(:,:,3) = sb_channel;
end
%--------------------------------------------------------------------------
function nrgb_img = normalized(rgb_img)
r_channel = rgb_img(:,:,1);
g_channel = rgb_img(:,:,2);
b_channel = rgb_img(:,:,3);

% nr_channel(r_channel ~= 0) = r_channel(r_channel ~= 0)/(r_channel(r_channel ~= 0) + g_channel(r_channel ~= 0) + b_channel(r_channel ~= 0));
% ng_channel(g_channel ~= 0) = g_channel(g_channel ~= 0)/(r_channel(g_channel ~= 0) + g_channel(g_channel ~= 0) + b_channel(g_channel ~= 0));
% nb_channel(b_channel ~= 0) = b_channel(b_channel ~= 0)/(r_channel(b_channel ~= 0) + g_channel(b_channel ~= 0) + b_channel(b_channel ~= 0));

if max(max(r_channel)) > 1.0 || max(max(g_channel)) > 1.0 || max(max(b_channel)) > 1.0
  r_channel = double(r_channel) / 255;
  g_channel = double(g_channel) / 255;
  b_channel = double(b_channel) / 255;
end

nrgb_img(:,:,1) = r_channel;
nrgb_img(:,:,2) = g_channel;
nrgb_img(:,:,3) = b_channel;
end